Codecov Report

❌ Patch coverage is 95.77465% with 3 lines in your changes missing coverage. Please review.

Re the test duplications, I think it could well make sense to do something similar to what GPUArrays does with its test suite: define compatible array and element type pairs and run the generically written testsuite on them. But that shouldn't hold this PR back, should we maybe open an issue for it?

WRT the tests, if qr and svd are supported, can higher level factorizations that use them (orthnull and buddies) also be tested with BigFloat?

In some tests, I had to change === to \approx

Does this mean that copies are happening where they shouldn't? What happens if you use == rather than ===?

I was under the impression that GenericLinearAlgebra also had a non-symmetric eigenvalue implementation but that might be wrong.

GenericLinearAlgebra supports the non-symmetric eigenvalue decomposition only when BigFloats aren't being used. For non-symmetric BigFloat matrices it throws an explicit error. In that regard, I would be fine with splitting this into two extensions.

Does this mean that copies are happening where they shouldn't? What happens if you use == rather than ===?

The == works, so indeed I'll check again whether I made some unintentional copies.

WRT the tests, if qr and svd are supported, can higher level factorizations that use them (orthnull and buddies) also be tested with BigFloat?

I'll take a look whether we can support them here. Probably the higher level factorizations should be fine, but things like the Schur decomposition or gen_eig might be less trivial.

I think I have now resolved all comments, apart from the comment from @Jutho about qr_Unblocked! and qr_Blocked!, and the comment from @lkdvos about the unnecessary allocations. I have opened an issue in both GenericLinearAlgebra.jl and GenericSchur.jl to support passing the output to these functions. Since we have an alternative (at least for the QR) with the new implemention of #90 that works generically and without extra allocations, I propose that we merge this PR in its current form. I can then also bring TensorKit.exp here, which was the ultimate goal. If the issues in the external packages are resolved, we can open a new PR to solve it here in the extensions as well.

If there are any additional comments, please let me know.

I'll take a look at these comments tomorrow. As for the allocations, I agree that the current implementation is suboptimal. My points was that this requires us indeed to intercept these allocations and change the tests (something also a little bit more annoying in the case of LQ, since we are using LQViaTransposedQR for this, and we have to intercept it internally), which is work that will just have to be reverted once the issues from the external packages are resolved.

I had a look at GenericSchur.jl, which I was much less familiar with than GenericLinearAlgebra.jl. That seems to reimplement a whole lot of LAPACK functionality related to Schur and eigenvalue decompositions, and so it might be possible to write the wrappers for schur_..., eig_ and eigh_ in a manner that is much closer in spirit to the LAPACK algorithms, and in particular writes directly to the given output arrays.

So I think it would definitely be good to split this PR into two parts corresponding to the two extensions.

GenericSchur.jl seems to provide both a QRIteration (which equals the current Francis name) and a DivideAndConquer implementation for the hermitian eigenvalue problem, and the standard QRIteration algorithm for the nonhermitian one (which for one or the other reason is not called as such in LAPACK). So we might want to think about the best naming strategy for these algorithms as well. I would definitely not call them GLA_... if they come from GenericSchur.jl, as that package might be worth while using in its own.

In case anyone is interested in reading more about the QR iteration method, Chapter 4 in https://people.inf.ethz.ch/arbenz/ewp/Lnotes/lsevp.pdf is really good.
My old implementation (which I still plan to include as well), is based on that.

I've now taken care of the 4 TODOs and the other comments.

I've renamed the algorithms to GLA_HouseholderQR, GLA_QRIteration, and GS_QRIteration to be as general and compatible with the other algorithms as possible. The folder is now split into two, namely genericlinearalgebra/ and genericschur/, i.e. one for each extension.

After a conversation with Jutho, we concluded that we think the best way forward is to resolve the two remaining points later. E.g.,

I had a look at GenericSchur.jl, which I was much less familiar with than GenericLinearAlgebra.jl. That seems to reimplement a whole lot of LAPACK functionality related to Schur and eigenvalue decompositions, and so it might be possible to write the wrappers for schur_..., eig_ and eigh_ in a manner that is much closer in spirit to the LAPACK algorithms, and in particular writes directly to the given output arrays.
So I think it would definitely be good to split this PR into two parts corresponding to the two extensions.

can be another PR where new algorithms are added, and

Considering the allocation of the output, do you know if we can pass the output to the implementations of these libraries? If not, that's definitely also okay, but in that case we should probably intercept allocating them in the first place, if they will be ignored. This can be achieved by overloading initialize_output(...) -> Nothing and altering the checks accordingly.

is now partly solved by using svd!, eigen!, eigvals!, and qr!. The problem will be solved completely once the other packages support passing the allocated output to their respective functions. I've opened an issue in both GenericSchur.jl (to which they have already answered) and GenericLinearAlgebra.jl, so hopefully this will be resolved relatively soon.

If anyone disagrees with this approach, or has some additional suggestions, please let me know.

I hope you don't mind, but I just went in and added the implementation of some of the comments that had already been discussed above and probably got overlooked, and I also simplified some of the code (based on suggestions that were already made above).

Some remarks for the implementations:

• svd has a full argument, which seems to be supported by GLA, and with the rmul!(one!(Q), Qprime) implementation we don't need the gramm schmidt procedure.
• I ended up not implementing/calling check_input, which I would argue is fair since these checks should already be handled by the library that we are wrapping. It could be straightforward to add this back if we feel this is necessary.

Would be good to get an additional review from @Jutho before merging, but this looks good to go for me.

Aside from the small comments above, I am a bit torn about is about registering these methods as the default_.... for StridedMatrix{BigFloat} and StridedMatrix{Complex{BigFloat}}.

On the one hand, they probably work for a larger set of types, such as Float16 eltype, or non-strided matrices.

On the other hand, especially with having in mind that we might have alternative Native_ algorithms for some of those, I am wondering if this is the best choice.

However, without these defaults in place, you would not be able to compute a matrix factorization with BigFloat entries without specifying the GLA_... algorithm explicitly, which is probably not what you want in generic code. If you explicitly load GenericLinearAlgebra, you probably want to use it.

I am however wondering what happens if you are using MatrixAlgebraKit and some other packages, some of which have GenericLinearAlgebra.jl as a dependency. Then also the package extension is loaded without explicit consent I assume?

I addressed the remainder of the comments since I was already busy anyways.

Considering the package dependency being loaded, I think your assumption is true, if anyone has this as a dependency somewhere, the code will be available and the default algorithms will be registered.
While I would in time also prefer to have our own native implementation as the default, I don't see too much harm in registering these default algorithms either.
We can still swap them afterwards?

What is interesting though is that GenericLinearAlgebra is definitely doing quite strong type piracy on the regular LinearAlgebra methods, and is actually simply hijacking them.
One alternative strategy therefore could be to, instead of extending GenericLinearAlgebra, simply extending LinearAlgebra itself, and using that as a fallback.
(IIRC, GenericLinearAlgebra.svd is actually just LinearAlgebra.svd imported).
There could be an argument about that being fair as a fallback too...

In any case, I don't think changing this in upcoming releases really matters too much, and even if this counts as a breaking change, it shouldn't really trickle too far downstream anyways.

I hope you don't mind, but I just went in and added the implementation of some of the comments that had already been discussed above and probably got overlooked, and I also simplified some of the code (based on suggestions that were already made above).

Not at all, thanks!

While I would in time also prefer to have our own native implementation as the default, I don't see too much harm in registering these default algorithms either.
We can still swap them afterwards?

I agree with Lukas here.

What is interesting though is that GenericLinearAlgebra is definitely doing quite strong type piracy on the regular LinearAlgebra methods, and is actually simply hijacking them. One alternative strategy therefore could be to, instead of extending GenericLinearAlgebra, simply extending LinearAlgebra itself, and using that as a fallback. (IIRC, GenericLinearAlgebra.svd is actually just LinearAlgebra.svd imported). There could be an argument about that being fair as a fallback too...

They are aware of this, and there is a PR to change it (JuliaLinearAlgebra/GenericLinearAlgebra.jl#147 ) but clearly that is not what people want, since indeed everyone wants the same call to just work using the fastest algorithm irrespective of the element type.

Ok, I think this is ready so I will merge.